

/*
program to calculate preferred habitat regressions

run_ph_regression
inputs: 
- first variable: shock

- mat_var: maturity remaining (years)
- ra_var: risk aversion high/low indicator
-- note: ra_var may be (semi-)cts: 0 in non-crisis, but cts in crisis
- long_var: term length long/short indicator
- security_id: variable with security identifier (eg cusip)
- auction_id: variable with auction identifier (eg datetime)

- tau_star: fixed maturity for rhs shocks (years)
- star_window: window for computing tau_star (years)

- cluster_type: T or NT for clustering on auction or two-way auction/security
-- optional: include bw(#) for HAC-robust SEs

- mat_window: size of maturity window to include in each regression (years)
- mat_long_window (optional): size of maturity window for maturities>20 
- mat_step: size of step forward while rolling over maturities (years)

- fname_stub (optional): save results in rolling data directory

outputs:
- results_low, results_high: matrices with regression output by maturity


*/

********************************************************************************
********************************************************************************
cap program drop run_ph_regression 
program define run_ph_regression

* parameter arguments
syntax varlist , ///
	mat_var(varname) ra_var(varname) long_var(varname) ///
	security_id(varname) auction_id(varname) ///
	tau_star(real) star_window(real) cluster_type(string) ///
	mat_window(real) mat_step(real) ///
	[mat_long_window(real 0) bw(real 0) fname_stub(string)]

***********************************
* first input variable: shock
local D_shock `varlist'
* create shock_star variables
* note: naming and referencing much easier using non-temp variables
cap drop D_star*

* tau_star maturity yield change
sort `auction_id' `mat_var' `security_id'
tempvar _shock_mean
bys `auction_id' (`mat_var' `security_id'): ///
	egen `_shock_mean' = mean(`D_shock') ///
	if `mat_var'>=`=(`tau_star' - 0.5*`star_window')' & ///
		`mat_var'<=`=(`tau_star' + 0.5*`star_window')'
bys `auction_id' (`mat_var' `security_id'): egen D_star = mean(`_shock_mean')
drop `_shock_mean'

* shock interactions (allowing for cts ra_var)
gen D_star_S0 = D_star*(`long_var'==0)*(`ra_var'==0)
gen D_star_S1 = D_star*(`long_var'==0)*(`ra_var')
gen D_star_L0 = D_star*(`long_var'==1)*(`ra_var'==0)
gen D_star_L1 = D_star*(`long_var'==1)*(`ra_var')

* create panel id: group security by short/long and low/high
tempvar panel_id
egen `panel_id' = group(`ra_var' `long_var' `security_id')

* maximum maturity included in regressions
summ `mat_var'
local mat_max `r(max)'
* for "rolling" windows: optionally increase for >20y maturities
if `mat_long_window' == 0 {
	local mat_long_window = `mat_window'
}
* final left point of maturity windows
local mat_window_final = `mat_max' - `mat_long_window'


***********************************
* variance-covariance options
* note: if bw specified, must use ivreghdfe (slower)
* xtset data (note: only matters when bw option specified)
tempvar T
sort `auction_id' `mat_var' `security_id'
egen `T' = group(`auction_id')
xtset `panel_id' `T'

if `bw'==0 {
	local reg_cmd reghdfe
	if "`cluster_type'"=="T" {
		* clustered on auction
		local vce_opts vce(cluster `T')
	}
	else if "`cluster_type'"=="NT" {
		* two-way clustered on auction/security
		local vce_opts vce(cluster `T' `panel_id')
	}
	else if "`cluster_type'"=="r" {
		* basic robust SEs
		local vce_opts vce(robust)
	}
}
else {
	local reg_cmd ivreghdfe
	if "`cluster_type'"=="T" {
		* clustered on auction
		local vce_opts cluster(`T') bw(`bw')
	}
	else if "`cluster_type'"=="NT" {
		* two-way clustered on auction/security
		local vce_opts cluster(`T' `panel_id') bw(`bw')
	}
	else if "`cluster_type'"=="r" {
		* basic robust SEs
		local vce_opts robust  bw(`bw')
	}
}
***********************************

* run regressions, rolling across maturities
local idx 0
local curr_min 0
while `curr_min' <= `mat_window_final' {
	n di "`curr_min'"
	local idx = `idx' + 1
	* include maturities between curr_min and curr_max
	local curr_min = (`idx'-1)*`mat_step'
	* optionally increased window size for long maturities
	if `curr_min' < 20 {
		local curr_max = `curr_min' + `mat_window'
	}
	else {
		local curr_max = `curr_min' + `mat_long_window'
	}
	* regression shock on shock star
	`reg_cmd' `D_shock' D_star_S0 D_star_S1 D_star_L0 D_star_L1 ///
		if `mat_var'>=`curr_min' & `mat_var'<`curr_max', ///
		a(`panel_id') `vce_opts'
	* save results to vector (dropping constant)
	mat b_vec = e(b)
	mat b_vec = b_vec[1, 1..4]
	* flatten variance-covariance matrix
	mat V_mat = e(V)
	mat V_mat = V_mat[1..4, 1..4]
	mata: st_matrix("V_vec", colshape( st_matrix("V_mat"), 16) )
	* get col/row names
	local V_rownames: rownames V_mat
	local V_colnames: colnames V_mat
	local V_vec_colnames ""
	forval r=1/4 {
		forval c=1/4 {
			local Vrn: word `r' of `V_rownames'
			local Vcn: word `c' of `V_colnames'
			local V_vec_colnames `"`V_vec_colnames' "V_`Vrn'_`Vcn'" "'
		}
	}
	matrix colnames V_vec = `V_vec_colnames'
	* additional stats
	* test equality of short/long coeffs
	test D_star_S0 = D_star_L0
	scalar p0 = r(p)
	test D_star_S1 = D_star_L1
	scalar p1 = r(p)
	mat reg_stats = (`curr_min', `curr_max', e(N), e(N_clust), ///
		e(df_r), e(r2_a_within), p0, p1)
	matrix colnames reg_stats = "maturity_start" "maturity_end" ///
		"N" "N_clust" "df_r" "r2_a_within" "eq_pval0" "eq_pval1"
	* combine
	mat res = (reg_stats[1, 1...], b_vec[1, 1...], V_vec[1, 1...])
	
	* start with row of all zeros
	if `idx'==1 {
		matrix results_mat = 0 * res
	}
	* append
	matrix results_mat = (results_mat \ res)
	
}
cap drop D_star*


* save results
if "`fname_stub'" != "" {
	preserve
	clear
	svmat results_mat, names(col)
	* save to dta file
	save ../data/rolling/phreg_`fname_stub'.dta, replace
	restore
}



end




